function [k,s,phi] = keffit1_1(c,b,a,dxx,nusig,sct,s,epsK,epsS )
% function [k,s,phi] = keffit(c,b,a,dxx,nusig,s,epsK,epsS )
%   c,b,a       = upper,lower,main diagonals of coeff matrix
%   dxx         = dxx volumes
%   nusig       = nu*sigmaf
%   s           = initial source guess
%   epsK,epsS   = convergence criteria for k and s
% uses the standard power iteration method to find k.

k       = 0.7;    % guess the initial keff
errK    = 1; 
errS    = 1; 
iter    = 0; 
itmx    = 1000;
sindx   = find(s);
n       = length(dxx);
numg    = length(nusig(1,:));
phi     = zeros(n,numg); 
kconv   = zeros(1,itmx); kconv(1,1) = k;
sconv   = zeros(n,itmx); sconv(:,1) = s;

while   ( ( errK>epsK ) || (errS>epsS) ) && iter < itmx
    s = s / sum(s); % is this a useful normalization?
    % call tridiagonal solver --------------------------
    phi(:,1)    = tridiag(b(:,1),c(:,1),a(:,1),s/k);
    scsrc       = zeros(n,1);
    if numg > 1                     % successive groups
        for g = 2:numg
            for gg = 1:(g-1)        
                for i = 1:n         % compute scattering source
                    scsrc(i,1) = scsrc(i,1) + sct(i,g,gg)*phi(i,gg);
                end
            end 
            % solve for phi_g
            phi(:,g)    = tridiag(c(:,g),b(:,g),a(:,g),s/k);
        end
    end
    sold = s; 
    kold = k;   
    s(:) = 0;       % reset fission source
    for g = 1:numg  % compute fission source
      s(:) = s(:) + nusig(:,g).*phi(:,g);
    end  
    % compute keff
    k = sum(s.*dxx)*kold/sum(sold(sindx).*dxx(sindx)); 
    
    % convergence stuff
%     kconv(1,iter+1)=k; 
%     sconv(:,iter+1)=s;
    errS = max( abs((s(sindx)-sold(sindx))./s(sindx)) );
    errK = abs( (k-kold)/k );
    iter = iter + 1;  % number of iterations
end

% some stuff for plotting convergence of k and s
% maxss = max(sconv);
% figure(2)
% subplot(2,1,1)
% plot(1:iter, kconv(1:iter),'kx')
% legend('k',1)
% subplot(2,1,2)
% plot(1:iter, maxss(1:iter),'r.' )
% legend('s',1)

disp(['*** final keff estimate: ', num2str(k),' ***'])
disp(['residual errors: errS = ',num2str(errS),' and errK = ',num2str(errK)])
disp(['in ',num2str(iter),' iterations'])
